home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / aminet / misc / sci / ephem_src_4_28.lha / precess.c < prev    next >
C/C++ Source or Header  |  1992-04-17  |  3KB  |  94 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. #define    DCOS(x)        cos(degrad(x))
  6. #define    DSIN(x)        sin(degrad(x))
  7. #define    DASIN(x)    raddeg(asin(x))
  8. #define    DATAN2(y,x)    raddeg(atan2((y),(x)))
  9.  
  10. /* corrects ra and dec, both in radians, for precession from epoch 1 to epoch 2.
  11.  * the epochs are given by their modified JDs, mjd1 and mjd2, respectively.
  12.  * N.B. ra and dec are modifed IN PLACE.
  13.  */
  14.  
  15. /*
  16.  * Copyright (c) 1990 by Craig Counterman. All rights reserved.
  17.  *
  18.  * This software may be redistributed freely, not sold.
  19.  * This copyright notice and disclaimer of warranty must remain
  20.  *    unchanged. 
  21.  *
  22.  * No representation is made about the suitability of this
  23.  * software for any purpose.  It is provided "as is" without express or
  24.  * implied warranty, to the extent permitted by applicable law.
  25.  *
  26.  * Rigorous precession. From Astronomical Ephemeris 1989, p. B18
  27.  */
  28.  
  29. precess (mjd1, mjd2, ra, dec)
  30. double mjd1, mjd2;    /* initial and final epoch modified JDs */
  31. double *ra, *dec;    /* ra/dec for mjd1 in, for mjd2 out */
  32. {
  33.     double zeta_A, z_A, theta_A;
  34.     double T;
  35.     double A, B, C;
  36.     double alpha, delta;
  37.     double alpha_in, delta_in;
  38.     double from_equinox, to_equinox;
  39.     double alpha2000, delta2000;
  40.  
  41.     mjd_year (mjd1, &from_equinox);
  42.     mjd_year (mjd2, &to_equinox);
  43.     alpha_in = raddeg(*ra);
  44.     delta_in = raddeg(*dec);
  45.  
  46.     /* From from_equinox to 2000.0 */
  47.     if (from_equinox != 2000.0) {
  48.         T = (from_equinox - 2000.0)/100.0;
  49.         zeta_A  = 0.6406161* T + 0.0000839* T*T + 0.0000050* T*T*T;
  50.         z_A     = 0.6406161* T + 0.0003041* T*T + 0.0000051* T*T*T;
  51.         theta_A = 0.5567530* T - 0.0001185* T*T + 0.0000116* T*T*T;
  52.  
  53.         A = DSIN(alpha_in - z_A) * DCOS(delta_in);
  54.         B = DCOS(alpha_in - z_A) * DCOS(theta_A) * DCOS(delta_in)
  55.           + DSIN(theta_A) * DSIN(delta_in);
  56.         C = -DCOS(alpha_in - z_A) * DSIN(theta_A) * DCOS(delta_in)
  57.           + DCOS(theta_A) * DSIN(delta_in);
  58.  
  59.         alpha2000 = DATAN2(A,B) - zeta_A;
  60.         range (&alpha2000, 360.0);
  61.         delta2000 = DASIN(C);
  62.     } else {
  63.         /* should get the same answer, but this could improve accruacy */
  64.         alpha2000 = alpha_in;
  65.         delta2000 = delta_in;
  66.     };
  67.  
  68.  
  69.     /* From 2000.0 to to_equinox */
  70.     if (to_equinox != 2000.0) {
  71.         T = (to_equinox - 2000.0)/100.0;
  72.         zeta_A  = 0.6406161* T + 0.0000839* T*T + 0.0000050* T*T*T;
  73.         z_A     = 0.6406161* T + 0.0003041* T*T + 0.0000051* T*T*T;
  74.         theta_A = 0.5567530* T - 0.0001185* T*T + 0.0000116* T*T*T;
  75.  
  76.         A = DSIN(alpha2000 + zeta_A) * DCOS(delta2000);
  77.         B = DCOS(alpha2000 + zeta_A) * DCOS(theta_A) * DCOS(delta2000)
  78.           - DSIN(theta_A) * DSIN(delta2000);
  79.         C = DCOS(alpha2000 + zeta_A) * DSIN(theta_A) * DCOS(delta2000)
  80.           + DCOS(theta_A) * DSIN(delta2000);
  81.  
  82.         alpha = DATAN2(A,B) + z_A;
  83.         range(&alpha, 360.0);
  84.         delta = DASIN(C);
  85.     } else {
  86.         /* should get the same answer, but this could improve accruacy */
  87.         alpha = alpha2000;
  88.         delta = delta2000;
  89.     };
  90.  
  91.     *ra = degrad(alpha);
  92.     *dec = degrad(delta);
  93. }
  94.